% =========================================================================
% modeFinderInitialize
% Run the MH simulator
% L. Melosi February 2018
%
% =========================================================================
c0=0;
c2=1e-5;

Nsim=5*1e5;
%% Uses modePosterior.m for the optimization
posStru.param.cal=posStru.param.cal;
posStru.steady.est=ssposest;

parStru.param=zeros(numpar, 1);
parStru.param(posStru.param.cal)=parStru.cal;

dataStru.data=Y;
dataStru.trainVec=trainvec;

filterStru.tauVec=addsol.tauVec;
filterStru.funcLikel=addsol.funcLikel;
filterStru.funcKfilter=addsol.funcKfilter;
filterStru.funcPost=@modePosterior;
filterStru.funcMin=@modePosteriorMin;

offset=1e-7;


%% posStru.param.estimated posStru.param.est
%% pos
%%
display('________________________________');
display('Begin Generating Starting Values');

npoints_full = 0;
% if priorStru.randomGridStart==0 && priorStru.randomPriorStart==0;
% use loopStru.parStartMat
disp('Using Starting values');
% loopStru.parStartMat=parstrvals;
% priorStru.numStartVals=size(loopStru.parStartMat,2);
%%

clear fcount npoints_full admissible lpostd two tone loop_logPost offset;
clear temp*;


%%
% =========================================================================
disp('Checking initial densities....');
[logPost,logLikel,logPrior]=...
    feval(filterStru.funcPost,loopStru.parMode(posStru.param.est,1),parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
if logPost < -1e15
    error('Likelihood too small! Parameter value is not admissible')
end

disp([sprintf('LogPosterior=%5.5e',logPost),sprintf('LogLikelihood=%5.5e',logLikel)]);
if abs( logPrior -( logPost - logLikel ) ) > 1e-3
    error('Log Posterior- Log Likelihood does not equal Log Posterior')
end

%Take the Hessian evaluated at the mode
HH=squeeze(loopStru.hessian(:,:,1));
paraold=loopStru.parMode(posStru.param.est,1);
%Draw an initial candidate
paraold=c0*chol(pinv(HH))'*randn(size(posStru.param.est,1),1)+paraold;
[logPost0,logLikel0,logPrior0]=...
    feval(filterStru.funcPost,paraold,parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
% logPost0=-logPost0;
   ii=0; paradraws=zeros(Nsim,size(posStru.param.est,1));
for i=1:Nsim
    %New draw
    flag=0;
    while flag==0
        paranew=c2*chol(pinv(HH))'*randn(size(posStru.param.est,1),1)+paraold;
        %Checking if the draws falls within bounds, if not, take another
        %one
        if sum(paranew<priorStru.ubnd(posStru.param.est))==size(posStru.param.est,1) && sum(paranew>priorStru.lbnd(posStru.param.est))==size(posStru.param.est,1)
            flag=1;
        end
    end
    [logPost1,logLikel1,logPrior1]=...
    feval(filterStru.funcPost,paranew,parStru,...
    posStru,model,dataStru,prpar,prss,filterStru,flags);
% logPost1=-logPost1;
    rr=exp(logPost1-logPost0);
    
    uu=rand(1,1);
    
    if uu<min(1,rr) %accept
        paraold=paranew;
        paradraws(i,:)=paranew';
        logPost0=logPost1;
        disp(sprintf('LogPosterior=%5.5e',logPost1))
        ii=ii+1;
    else  %Discard the draw
        paradraws(i,:)=paraold';
    end
    disp(sprintf('Number of Draws Done So Far'))
    i
    disp(sprintf('Acceptance Ratio'))
    ii/i*100
    
end
save PostDraws/postDraws paradraws